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ABSTRACT 

We present the results of three-dimensional hydrodynamical simulations of the final 
stages of inspiral in a black hole-neutron star binary, when the separation is compara- 
ble to the stellar radius. We use a Newtonian Smooth Particle Hydrodynamics (SPH) 
code to model the evolution of the system, and take the neutron star to be a polytrope 
with a soft (adiabatic index F = 2 and F = 5/3) equation of state and the black hole 
to be a Newtonian point mass. The only non-Newtonian effect we include is a gravi- 
tational radiation back reaction force, computed in the quadrupole approximation for 
point masses. We use irrotational binaries as initial conditions for our dynamical sim- 
ulations, which arc begun when the system is on the verge of initiating mass transfer 
and followed for approximately 23 ms. For all the cases studied we find that the star 
is disrupted on a dynamical time-scale, and forms a massive {Mdisc ~ 0.2 Mq) accre- 
tion torus around the spinning (Kerr) black hole. The rotation axis is clear of baryons 
(less than 1Q~^Mq within 10°) to an extent that would not preclude the formation 
of a relativistic fireball capable of powering a cosmological gamma ray burst. Some 
mass (the specific amount is sensitive to the stiffness of the equation of state) may 
be dynamically ejected from the system during the coalescence and could undergo r- 
process nucleosynthesis. We calculate the waveforms, luminosities and energy spectra 
of the gravitational radiation signal and show how they reflect the global outcome of 
the coalescence process. 

Key words: binaries: close — gamma rays: bursts — gravitational waves — hydro- 
dynamics — stars: neutron 



1 INTRODUCTION AND MOTIVATION 

In binary systems, the emission of gravitational waves and 
accompanying loss of angular momentum will lead to a de- 
crease in the orbital separation, and coalescence will occur if 
the decay time is small enough (less than the Hubble time). 
For binaries made of neutron stars, PSR 1913-f 16 being the 
most famous example, this consequence of general relativ- 
ity has been observed indirectly (Hulse & Taylor 1975) (see 
also Wolszczan (1991) for the case of PSR 1534-1-12), and 
the change in orbital period matches the theoretical predi- 
cion to very high accuracy (Taylor et al. 1992; Stairs et al. 
1998). Given their present-day orbital periods (on the order 
of 10 hours), these systems will eventually merge. The final 
stages of the coalescence present an opportunity to study 
the equation of state at very high densities (the system 
is in effect a giant accelerator), and will undoubtedly pro- 
duce a strong electromagnetic and gravitational wave signal 



containing some of this information. No black hole-neutron 
star binary systems are known yet, but population synthe- 
sis studies (Lattimer & Schramm 1976; Narayan, Piran & 
Shemi 1991; Tutukov & Yungelson 1993; Lipunov, Postnov 
& Prokhorov 1997; Bethe & Brown 1998; Portegies Zwart 
& Yungelson 1998; Belczyhski & Bulik 1999; Kalogera et 
al. 2001; Kalogera & Belczyhski 2001) over the past 25 
years lead one to believe that their rate is comparable to 
that of double neutron star binaries, and is on the order of 
10"'' — 10~^ per year per galaxy. 

Solving this problem completely clearly requires de- 
tailed hydrodynamic modeling in three dimensions, radi- 
ation transport, a realistic equation of state, and general 
relativity. As such, it must be approached in stages, with 
successive approximations depending on the aspect of the 
general problem one wishes, and is able, to solve. 

Compact binaries are expected to be sources of grav- 
itational radiation observable by detectors such as LIGO 
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(Abramovici et al. 1992) and VIRGO (Bradaschia ct al. 
1990) as the inspiral occurs. The signal can be approxi- 
mated as that of point-masses and calculated accurately us- 
ing post-Newtonian expansions when the separation is large, 
compared with the stellar radius (Kidder, Will & Wiseman 
1992; Cutler et al. 1994; Blanchet et al. 1995). When the dis- 
tance between the stars becomes comparable to their radii, 
hydrodynamical modeling becomes essential. The theoret- 
ical study of the tidal disruption of a neutron star by a 
black hole was addressed many years ago (Wlicclcr 1971; 
Lattimer & Schramm 1974; Lattimer & Schramm 1976), 
and numerical hydrodynamical studies of binary neutron 
star coalescence were begun somewhat more recently, ini- 
tially by Oohara & Nakamura (1989; 1990; 1992) and Naka- 
mura & Oohara (1989; 1991). The work of Chandrasekhar 
(1969) on incompressible ellipsoids was generalized to the 
compressible case in the Newtonian regime by Lai, Rasio & 
Shapiro (1993b, hereafter LRSb), using a polytropic equa- 
tion of state, who showed that tidal effects alone could pro- 
duce a de-stabilization of the orbit in certain situations (Lai, 
Rasio & Shapiro 1993a). Rasio & Shapiro (1992; 1994; 1995, 
hereafter RS92, RS94, RS95 respectively) then performed 
a series of dynamical simulations to study the coalescence 
of two neutron stars, using Smooth Particle Hydrodynam- 
ics (SPH), while Zhuge, CentrcUa & McMillan (1994; 1996) 
focused on the gravitational waves spectrum. Both of these 
groups used a polytropic equation of state throughout their 
analysis. The thermodynamical details of the process were 
studied by Davies et al. (1994), Ruffert, Janka & Schafer 
(1996), Ruffert et al. (1997) and Rosswog et al. (1999; 2000), 
by using the equation of state of Lattimer & Swesty (1991). 
This work was all done using a Newtonian or modified New- 
tonian approach (by including gravitational radiation reac- 
tion in different ways in the calculations), and we note that 
the thermodynamic details are of little importance for the 
emission of gravitational waves, since it is conccrnccd pri- 
marily with the motion of bulk matter at high densities. 
More recently, there have been advances iu making post- 
Newtonian calculations of initial conditions (Lombardi, Ra- 
sio & Shapiro 1997) and mergers (Faber & Rasio 2000; 
Faber, Rasio & Manor 2000; Ayal et al. 2001), and also in 
including general relativity (Wilson, Mathews & Marronetti 
1996; Baumgartc et al. 1997; Oohara & Nakamura 1999; 
Shibata 1999; Shibata k Uryu 2000; Uryii & Eriguchi 2000; 
Usui, UryiJ & Eriguchi 2000; Gourgoulhon et al. 2001). 

The gamma ray hurts (GRBs) are now believed to be 
at cosmological distances (Meegan et al. 1992), after the 
discovery of optical afterglows (Meszaros & Rees 1997a) in 
the last few years that have established their redshifts (Met- 
zger et al. 1997; Djorgovski et al. 1998; Kulkarni et al. 1998; 
Kulkarni et al. 1999). Reviews have been given by Fishman 
& Meegan (1995) and van Paradijs, Kouveliotou & Wijers 
(2000). Observations have shown that (i) there is a bimodal- 
ity in burst durations (Kouveliotou et al. 1995), with classes 
of short {tburst — 0.5 s) and long {tturst — 40 s) bursts; 
(ii) the energy release if one assumes isotropy is on the or- 
der of 10^^ — 10^^ erg; (iii) many bursts show variability 
on very short (ms) time-scales; (iv) at least some bursts 
arc beamed, implying a lower energy release than isotropic 
emission would lead one to believe (see e.g. Harrison et 
al. (1999); Stanek et al. (1999); Frail et al. (2001)). The pre- 
ferred model for GRBs involves the expansion of a relativis- 



2.35 




2.00 - / 
1.95 - / //' 
1.90 - ' '/ 
1.85 / 

1.80 L- ^ ' ' ' ' ' ' 

2.6 2.8 3 3.2 3.4 3.6 3.8 4 
r 




2.8 - / y-y 

2.7 ^ ' ' ' < ' ' ' ' ' 

2.9 3 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 4 
r 

tic fireball (Rees & Meszaros 1992; Meszaros & Rees 1993) 
which would produce the gamma rays through relativistic 
shocks and subsequent synchrotron radiation. The fireball 
would presumably originate from a central engine capable 
of accomodating the observational requirements mentioned 
above. A variety of progenitors involving compact objects 
have been suggested, see e.g. Fryer, Woosley & Hartmann 
(1999a). Many of them invoke an accretion torus around a 
black hole, originating from a double neutron star coales- 
cence (Goodman 1986; Paczyriski 1986; Eichler et al. 1989; 
Narayan, Paczyhski & Piran 1992), where the central object 
would presumably collapse to a black hole, the merger of a 
neutron star, white dwarf or Helium core with a black hole 
(Paczyriski 1991; Fryer et al. 1999b; Zhang & Fryer 2001), or 
a "failed supernova" or coUapsar (Woosley 1993; MacFadyen 
& Woosley 1999), where a massive star collapses but pro- 
duces a black hole instead of a neutron star at its center. 
These systems would produce a GRB by tapping the bind- 
ing energy of the disc through neutrino emission (Goodman 
1987; Jarozynski 1993; Jarozynski 1996; Witt et al. 1994; 
Thompson 1994; Mochkovitch et al. 1993; Mochkovitch et 
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Figure 1. Total angular momentum as a function of binary sepa- 
ration for irrotational Roche— Riemann binaries using the method 
of LRSb for r = 2 (solid lines) and T = 5/3 (dashed lines) for (a) 
q=0.5; (b) q=0.31; (c) q=0.2. The dotted lines correspond to the 

solution computed for point masses in Kcplcrian orbits. The thick 
vertical lines mark the initial separations used for the dynamical 
runs. 

al. 1995; Popham, Woosley & Fryer 1999; Ruffert & Janka 
1999), or the rotational energy of the black hole, through the 
Blandford & Znajck (1977) mechanism, producing so-called 
Poynting jets (Meszaros & Rees 1997b; Meszaros, Rees & 
Wijers 1999; Lee, Wijers & Brown 2000). Another class of 
models also involves neutron stars, but would power the 
GRB through the catastrophic release of rotational energy 
via intense magnetic fields (Usov 1992; Kluzniak & Ruder- 
man 1998; Spruit 1999; Ruderman, Tao & Kluzniak 2000), 
or even through intense neutrino emission in a neutron star 
binary before the coalescence, because of tidal heating and 
compression (Salmonson, Wilson & Mathews 2001). 

The ejection of neutron star matter to the interstellar 
medium during a dynamical coalescence might contribute 
to the abundances of heavy elements (Lattimer & Schramm 
1974; Lattimer & Schramm 1976; Symbalisty & Schramm 
1989; Eichler et al. 1989), in addition to the amounts ex- 
pected from supernova explosions (Meyer & Brown 1997; 
Preiburghaus et al. 1999a). This question has been addressed 
in the numerical calculations of double neutron star merg- 
ers of Rosswog et al. (1999; 2000) and by Preiburghaus et 
al. (1999b). If the rates of black hole-neutron star mergers 
are comparable, it is possible that these systems might also 
contribute in the same way to the galactic abundances. 

Our work on merging black hole-neutron star bina- 
ries began with low-resolution simulations (Lee & Kluzniak 
1995) that used a stiff polytropic equation of state. The re- 
sults initially led us to believe that, if proved true, these 
models were not likely to produce cosmological gamma ray 
bursts, because of excessive baryon contamination. We grad- 
ually increased our immerical resolution, using essentially 
Newtonian physics (except for our treament of gravita- 
tional radiation reaction, see below, § 2), and treated tidally 
locked binaries with stiff and soft equations of state (Lee & 



Kluzniak 1999a; Lee & Kluzniak 1999b, hereafter papers I 
& II), and irrotational binaries with a stiff equation of state 
(Lee 2000, hereafter paper III), always using a polytrope to 
model the initial neutron star. It became apparent early on 
that our initial suspicions were unfounded, and that indeed 
these systems were promising candidates for the central en- 
gines of GRBs (Lee & Kluzniak 1997; Kluzniak & Lee 1998). 
We also found that for very stiff equations of state (see pa- 
pers I and III) the neutron star could avoid immediate tidal 
disruption, and that this would be reflected in the gravi- 
tational wave signal. Additionally, a substantial amount of 
mass could be ejected to the interstellar medium, and po- 
tentially undergo r-process nucleosynthesis, thus contribut- 
ing to the abundances of heavy elements. Recently, Janka 
et al. (1999) used the same formalism that Ruffert et al. 
(1996; 1997) had employed for binary neutron star coales- 
cence studies to simulate the merger of a black hole with 
a neutron star. Their calculations have revealed the same 
qualitative aspects of the process which we have found, with 
differences due mainly to the different formalism used for 
gravitational radiation reaction and their use of a different 
equation of state (Lattimer & Swesty 1991). 

This paper is the last in the series that has used the ap- 
proach briefly described above (and detailed below in § 2), 
having mapped the parameter space we intended to explore 
by varying the stiffness of the equation of state, the initial 
mass ratio in the binary and the distribution of angular mo- 
mentum in the system (using tidally locked and irrotational 
binaries as initial conditions). A short exposition on the nu- 
merical method and initial conditions is given in § 2 and 
§ 3 (for details concerning the implementation we refer the 
reader to the longer corresponding sections in paper III and 
the appendix in paper I), followed by our results in § 4. The 
effect different choices of initial conditions can have on the 
dynamical coalescence is presented in § 5, and a summary 
and discussion is given in § 6. 



2 NUMERICAL METHOD 

For the calculations presented in this work, we have used the 
method known as Smooth Particle Hydrodynamics (SPH) 
(see Monaghan 1992 for a review and Lee 1998 for a de- 
scription of our own code). The code is the same one that 
was used for our previous simulations of irrotational black 
hole-neutron star binaries (paper III). Here we will not dis- 
cuss the code in detail, but limit the presentation to a few 
basic points. 

As before, the black hole is modeled as a Newtonian 
point mass of mass Mbh with an absorbing boundary at the 
Schwarzschild radius rsch = 2GMbh/c^ ■ Any SPH particle 
that crosses this boundary is absorbed by the black hole, 

whose mass and momentum are adjusted so as to ensure 
conservation of total mass and total linear momentum in 
the system. 

The neutron star is modeled as a polytrope with a soft 
equation of state, so that the pressure is given by P = Kp^ 
with r and K being constants (see paper I). Unless oth- 
erwise noted, we measure mass and distance in units of the 
mass and radius of the unperturbed (spherical) neutron star 
(13.4 km and 1.4 M© respectively), so that the units of time. 
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density and velocity are 




For the dynamical simulations presented here, we have used 
N ~ 80, 000 SPH particles to model the neutron star. The 
initial (spherical) polytropo is constructed by placing the 
SPH particles on a uniform three-dimensional grid with par- 
ticle masses proportional to the Lane-Emden density. This 
ensures that the spatial resolution is approximately uniform 
throughout the fluid. This isolated star is then allowed to 
relax for a period of thirty time units (as defined above) 
by including a damping term linear in the velocities in the 
equations of motion. The specific entropies of the particles 
are kept constant during this procedure (i.e. K is constant 
in the equation of state P — Kp^). 

To perform a dynamical run, the black hole and ev- 
ery SPH particle are given the velocity as determined from 
the corresponding initial condition (see below) in an inertial 
frame, with the origin of coordinates at the centre of mass of 
the system. Each SPH particle is assigned a specific internal 
energy Ui = Kp^^~^^ /{T — 1) and the equation of state is 
changed to that of an ideal gas, P = (T — l)pu. The specific 
internal energy is then evolved according to the first law of 
thermodynamics, taking into account the contribution from 
artificial viscosity (see below) . Wc vary the initial mass ratio 
q = Mj^s/Mbh in the binary by adjusting the mass of the 
black hole only. 

Artificial viscosity is used in SPH to handle the pres- 
ence of shocks and avoid particle interpenetration. As in pa- 
per HI, we use the form of Balsara (1995), which vanishes in 
regions of large vorticity but retains the ability to deal with 
the presence of shocks (in regions of strong compression). 

We include gravitational radiation reaction in the 
quadrupole approximation for point masses (Landau & Lif- 
shitz 1975), with the same implementation as described in 
paper HI. Namely, we apply a back reaction force to the 
black hole and the self bound core of the neutron star, 
treating the latter as a point mass. The corresponding 
terms in the equations of motion axe switched off once the 
star is tidally disrupted, when the core mass drops below 
0.14 Mq. We continuously compute the radiation reaction 
time-scale t'l^ = 256G^MBHMcorc(MBH + M^ore)/{5r*c^) 
and an estimate of the current orbital period torb = 
27r/ G(Mbh + Mcoie)/r^, where r is the separation be- 
tween the black hole and the centre of meiss of the core. For 
the typical separations and masses in the black hole-core 
binary, by the time the core mass has dropped to 0.14 Mq, 
the radiation reaction time-scale is much longer (by at least 
an order of magnitude) than the current orbital period. 



3 INITIAL CONDITIONS 

Exactly as for the results presented in paper HI, we have 
used irrotational binaries for the dynamical runs shown here. 
This amounts to considering the stars have no spin in an 
external, inertial frame of reference. This initial condition 
is more realistic than that of a tidally locked binary, be- 
cause the viscosity inside neutron stars is not large enough 
to maintain synchronization during the inspiral phase (Bild- 
sten & Cutler 1992; Kochanek 1992). Essentially, the stars 
will coalesce with whatever spin angular momenta they have 
when the binary separation is large. Setting up accurate and 
self-consistent initial conditions in this case is not a triv- 
ial matter, and we use the same approximation as before. 
Namely, we apply the energy variational method of LRSb 
and take the neutron star to be a compressible tri-axial 
Roche-Riemann ellipsoid (see § 8 in LRSb). 

We build an initial condition by first constructing a 
spherical star of given radius and mass, as described in § 2. 
We then use the method of LRSb to calculate the orbital 
angular velocity of the binary and the semimajor axes of 
the Roche-Riemann ellipsoid for the appropriate choice of 
adiabatic index, initial mass ratio and binary separation (see 
Table 1). The semimajor axes of the fluid configuration can 
also be calculated from the SPH numerical solution using 

/ 57« 

V KnMNS 

where 

i 

The stifi^ncss of the equation of state enters these equa- 
tions through the parameter (/t„ = 0.653 for F = 2 and 
Kn = 0.511 for F = 5/3). The first and second semimajor 
axes of the tri axial ellipsoid lie in the orbital plane, with 
the first one along the line joining the two binary compo- 
nents. The third a^cis is oriented perpendicular to this plane 
(along the axis of rotation). The position of each SPH par- 
ticle is then re-scaled (independently along each coordinate 
axis) so that the new fluid configuration has the appropri- 
ate semimajor axes. This ellipsoid is then used as an initial 
condition for the corresponding dynamical run. The initial 
velocity is given by the orbital angular velocity (for the az- 
imuthal component) plus the radial velocity corresponding 
to point-mass inspiral. The variation in total angular mo- 
mentum as a function of binary separation for irrotational 
Roche-Riemann binaries (with various mass ratios and adi- 
abatic indices in the equation of state) is shown in Figure 1, 
as calculated using the method of LRSb. The curves show 
a turning point as the separation is decreased, indicating 
the presence of a dynamical instability in the system. Two 
distinctions are necessary at this point. First, the ellipsoidal 
approximation becomes less accurate as the separation is 
decreased. This applies regardless of the value of the adia- 
batic index, but is much more serious for stiff equations of 
state, because the tidal effects are more pronounced. Sec- 
ond, the adiabatic index does determine if the Roche limit 
(when mass transfer through overflow of the lobe occurs) is 
reached before or after the dynamical instability. For stiff 
equations of state (such as the ones shown in papers I & 
III), the instability can be reached at or before the Roche 



© 2001 RAS, MNRAS 000, 1-21 



Black hole-neutron star coalescence 5 



limit. However, for a more compressible case (see paper II) 
the inverse occurs, and it is the mass transfer process itself 
(which is unstable) that is responsible for the subsequent 
evolution of the system. 

Wc have chosen the values of the initial separation for 
our dynamical runs r, to be slightly above the turning point 
(see Table 1). The ellipsoidal approximation is then still rear 
sonable, and our equilibrium configurations have not yet 
reached the point where the neutron star will overflow its 
Roche lobe. When the dynamical simulation is initiated, 
the separation will decrease due to the emission of gravi- 
tational waves, and mass transfer will start promptly. The 
construction of full equilibrium initial conditions at the point 
of Roche lobe overflow is a problem that was addressed by 
Uryti & Eriguchi (1999). 

The initial separations we have chosen are similar to 
what we have presented before for the case of tidally locked 
(papers I & II) and irrotational (paper III) black hole- 
neutron star binaries. We also include in Table 1 the initial 
parameters for two runs (C31S and D31S) that used initially 
spherical neutron stars for the dynamical calculations, with 
a Keplerian orbital angular velocity (as for runs A31S and 
B31S in paper III). We have performed these runs to gauge 
the effect non-equilibrium initial conditions will have on the 
evolution of the system, and show the results in § 5. 

4 RESULTS 

4.1 Morphology of the mergers 

For every dynamical run, the decrease in binary separation 
leads to Roche lobe overflow on an orbital time-scale. A 
stream of gas forms at the inner Lagrange point, transfer- 
ring matter from the neutron star to the black hole. At the 
same time, the star is tidally stretched and extends away 
from the black hole through the external Lagrange point. 
We show in Figures 2 and 3 density contours in the orbital 
plane at various times during the simulation for runs C31 
and D31. For aU other runs (C51, C31S, C20, D51, D31S 
and D20) the plots are qualitatively similar. As the accre- 
tion stream winds around the black hole, it collides with it- 
self and forms a torus, while the gas thrown out through the 
outer Lagrange point forms a long tidal tail. For a the less 
compressible case (F = 2, run C31), the torus, as well as the 
tidal tail, are thinner, as one should expect. The accretion 
torus that is formed around the black hole is not initially 
azimuthally symmetric, but shows a double ring structure, 
particularly for F = 2 (see panels (d)-(h) in Figure 2). This 
appears as the gas that passes through periastron near the 
black hole overshoots the circular orbit that would corre- 
spond to the specific angular momentum it contains, form- 
ing an outer ring (see panels (c)-(d) in Figure 2). It then 
falls back towards the black hole and encounters the rear 
of the accretion stream. The subsequent collision tends to 
circularize the orbit of the fluid, and also pushes it to the 
inner ring, closer to the black hole (panels (d)-(e)). The 
structure of the outer ring rotates slowly counterclockwise 
(with the initial orbital motion) as the bulk of the tidally 
disrupted star (which produces the accretion stream) con- 
tinues orbital motion in the same direction, on the opposite 
side of the black hole. At late times, the density contrast be- 
tween the rings drops (see below. Figure 7(a) and (c), and 



Figure 10), but nevertheless a hump remains in the accre- 
tion disc, as there is still a visible stream feeding it from the 
oposite side. This structure was clearly seen for F = 2.5 in 
the results presented in paper III, and for the same reasons. 
It is present as well for F = 5/3 (see Figure 3), although the 
distinction between having two rings and a hump is not as 
marked, even as the disc is forming (panels (c)-(e)). This is 
due to the higher compressibility of the material, and hence 
its tendency to expand at low densities. At late times the 
disc is much more azimuthally symmetric than for F = 2 
(see below Figure 7(b) and (d)). 

Our implementation of gravitational radiation reaction 
is valid only for circular orbits. Thus wo monitor the ec- 
centricity e of the orbital motion during the coalescence, to 
ensure that it remains small before gravitatational radia- 
tion reaction is switched off. We compute an estimate for e 
by assuming that it is that of a binary system with masses 
Mbh and Mcore, given by e = y/l + 2EJ'2/G'^^iM?orcM^^, 
where fj, = McoreAfBH/(Mcore -l- Mbh), and E and J are the 
mechanical energy and angular momentum of the orbital 
motion. During the initial phase, e ~ 0.05, and close to the 
instant of minimum binary separation (at t ~ 20 for most 
runs) e < 0.1. By the time radiation reaction is switched 
off (at t ~ 30, see Table 1), the eccentricity has increased 
somewhat, to e ~ 0.2. At this stage the mass ratio (between 
the core and the black hole) has dropped enough so that the 
effects of including radiation reaction are very small. 

The separation between the centre of mass of the core 
and the black hole is shown in Figure 4 for all runs. Initially, 
it decreases at a rate consistent with that of a point-mass 
binary, and subsequently does so at an even faster rate, due 
to hydrodynamical effects. This is particularly important for 
high mass ratios {q — 0.5 and q = 0.31). For q = 0.2, the 
deviation is smallest and almost negligible, until t ~ 13—15, 
depending on the value of F. Once the initial mass transfer 
episode is under way, the separation reaches a minimum and 
then increases, as the core of the star is stretched and moves 
to a greater separation. Qualitatively, the evolution resem- 
bles that of a stiff equation of state (paper III) except that 
now the point-mass approximation for the orbital decay is 
valid for a longer time (at smaller separations for a given 
value of q, compare for example the case with q = 0.2 in the 
two panels in Figure 4). This is simply due to the fact that 
the stars are better approximated by point masses as the 
adiabatic index is decreased. The point at which the separar 
tion is at a minimum coincides with the maximum accretion 
rate (sec below). After this initial periastron passage, the 
star is completely disrupted, and in every case the flnal con- 
figuration consists of a massive accretion disc around the 
black hole. Gravitational radiation reaction is switched off 
at t ~ 30 for all runs (see Table 1) when the core mass drops 
below 0.14 Mq. 

As mentioned above, the material which moves away 
from the black hole through the outer Lagrange point forms 
a long one-armed spiral in the system. This structure is 
usually formed during a dynamical coalescence (it is a two- 
armed spiral in the case of neutron star mergers, with each 
star producing one arm, sec e.g. RS94). The main difference 
between the runs presented here and the case of low com- 
pressibility studied in paper III is that no clumps are formed, 
and the distribution of matter remains smooth along the 
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Table 1. Basic parameters for each run. The table lists for eax;h run (labeled) the initial mass ratio, the adiabatic index used, the 
initial orbital separation, the axis ratios for the tri-axial ellipsoid used as an initial condition, the initial orbital angular velocity of the 
binary, the time at which gravitational radiation rea<;tion is switched off in the simulation, the time at which the simulation was stopped, 
and the initial number of particles. The runs labeled C31S and D31S used an initially spherical neutron star (otherwise irrotational 
Roche-Riemann ellipsoids were used, see text for details). 



Run 
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02/01 


as/ai 
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^rad 


*/ 


N 


C50 


0.50 


2.0 


3.25 


0.842 


0.857 


0.29753 


35.12 


200.0 


81,608 


C31 


0.31 


2.0 


3.70 


0.828 


0.844 


0.29042 


35.80 


200.0 


81,608 


C31S 


0.31 


2.0 


3.70 


1.000 


1.000 


0.28881 


34.76 


200.0 


81,608 


C20 


0.20 


2.0 


4.15 


0.808 


0.826 


0.29117 


30.12 


200.0 


81,608 


D50 


0.50 


5/3 


3.25 


0.904 


0.911 


0.29644 


37.03 


200.0 


82,136 


D31 


0.31 


5/3 


3.60 


0.884 


0.892 


0.30182 


30.40 


200.0 


82,136 


D31S 


0.31 


5/3 


3.60 


1.000 


1.000 


0.30094 


29.44 


200.0 


82,136 


D20 


0.20 


5/3 


4.15 


0.884 


0.892 


0.29037 


26.96 


200.0 


82,136 



length of the tail (as observed also in paper II, RS92, RS95 
for soft equations of state). 

Wc show in Figure 5 the accretion rates onto the black 
hole for runs C50, C31, C20, D50, D31 and D20. The 
maximum rates are reached during the initial episode of 
Roche lobe overflow {Mmax ~ 0.06 — 0.09, equivalent to 
0.7 — 1.1 Mq ms~^). They correspond mainly to matter 
that is directly accreted by the hole from the mass transfer 
stream coming from the neutron star. As the accretion disc 
is formed around the black hole, M gradually decreases, al- 
though there are small oscillations. One can see in the curves 
that there are secondary maxima in M at f < 60 for all runs. 
This is due to the circularization process of the orbits in the 
disc. The streams of matter coming from the neutron star 
collide with themselves, and some matter (along the inner 
edge of the stream) falls onto the black hole with greater 
ease, giving rise to the quasi-periodic oscillations in the ac- 
cretion rate. This only occurs two or three times at most, 
and at late times {t > 100) the accretion rate is decreas- 
ing monotonically, showing the circularization of the orbits. 
These peaks are present in the runs shown in paper III, but 
the lower resolution used there makes it harder to appreciate 
them. 

The peak accretion rates shown in Table 2 are substan- 
tially higher when F = 5/3, by up to a factor of ~ 1.4 for 
q — 0.2. This is one of many effects of the mass radius re- 
lationship that are observed during dynamical coalescences. 
For polytropes, R oc m'^"^/^^^"*', so for F = 2 the stellar 
radius is unaffected by mass loss (or accretion), while for 
F = 5/3, R oc M~^^^ and thus the star will expand upon 
overflowing its Roche lobe and losing mass. The decrease 
in separation due to energy losses to gravitational waves is 
what initially brings the star to the point of Roche lobe over- 
flow. The star's reaction when this happens then depends on 
the compressibility. When F = 2, the ofi^ect of tidal forces 
and gravitational wave back reaction are enough to com- 
pletely disrupt the star. For F = 5/3, there is an added, 
runaway effect, because since the star expands upon losing 
mass, it further overflows its Roche lobe (and thus produces 
higher accretion rates). This alone can de-stabilize the or- 
bit and induce coalescence, as observed for the case of a 
tidally locked system with F = 5/3, without gravitational 
wave back reaction (this was done initially as a test and re- 
ported in paper II). This also explains why the total disc 
mass is lower for a given initial mass ratio for F = 5/3 (also 
given in Table 2, column 2). 

The total angular momentum in the system decreases 



for two reasons. First, there is a decay due to the emission 
of gravitational waves (seen at early times before substan- 
tial mass transfer has taken place), and second, much of 
the angular momentum of the accreted matter is lost to the 
spin of the black hole. As stated above, when accretion oc- 
curs we update the mass and momentum of the black hole 
so as to ensure conservation of total mass and linear mo- 
mentum. Conservation of angular momentum then allows 
us to estimate the degree to which the black hole is spun up 
as a result of accretion, and we calculate its Kerr parameter 
a — Jbhc/GMIpj at the end of the simulation (we take a = 
at t = 0). This is shown in column 8 of Table 2. We note 
that our previously published results for a in papers II and 
III contained an error which we have now corrected. An ex- 
planation and the correct values are given in the appendix. 
For a fixed adiabatic index, the black hole is spun up to a 
greater degree (up to almost one half the maximum rotation 
rate) at higher mass ratios, simply because it is less massive. 
The higher rotation rates seen at lower F reflect the corre- 
sponding higher accretion rates (see above) and the fact that 
the total accreted mass is greater. 

The various energies in the system are shown in Fig- 
ure 6 for runs C31 and D31. The changes in mechanical 
energy seen at early times are due to the back reaction of 
gravitational waves, while the initial episode of mass trans- 
fer in the initial stages of the coalescence is evident in the 
large changes that occur at t « 15. At later times, the vari- 
ation is small and the curves show a monotonical decay as 
the accretion discs become more azimuthally symmetric. 

The core of the neutron star moves away from the black 
hole for the same reasons as outlined in paper III for the 
case of a stiff equation of state. In the case of conservative 
mass transfer (where the total mass and orbital angular mo- 
mentum J are conserved), if the donor is the less massive 
component, the binary separation will increase. The system 
is not strictly conservative in this case, but the global re- 
sponse is the same. The specific angular momentum in the 
core increases as mass transfer proceeds, and this makes the 
separation increase. The mass-radius relationship outlined 
above makes it impossible for the system to survive as a sta- 
ble binary, as was the case for F = 3 (paper III). As soon as 
the star overflows its Roche lobe, catastrophic mass transfer 
ensues, completely disrupting the star. The two dominant 
effects as far as the orbital evolution is concerned are the 
gravitational wave emission (and the accompanying loss of 
angular momentum) and mass transfer. In the cases shown 
here, (where the neutron star expands or maintains a con- 
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Figure 2. Density contours in tlie orbital plane during the dynamical simulation of the black hole-neutron star binary with initial mass 
ratio q = 0.31 and F = 2 (run C31). The orbital rotation is counterclockwise. All contours are logarithmic and equally spaced every 0.25 
dex. Bold contours are plotted at logp = —4, —3, —2, —1 (if present) in the units defined in eq. 2. The thick black arcs bound the matter 
that forms the core (see § 2). The time for each frame is given in the units defined in eq. 1. 



stant radius upon losing mass), both effects lead to complete 
tidal disruption of the star on an orbital time-scale. For a 
stiff equation of state they tended to drive the system in 
opposite directions, with angular momentum losses making 
the separation decrease while mass transfer increased it. The 
outcome in that case was episodic mass transfer, and the 
frequent formation of accretion discs when the star was dis- 
rupted. Thus, in the present case as well, stable mass transfer 
from the neutron star is impossible, and the final configu- 
ration consists of a massive accretion disc around the black 
hole. 



4.2 Accretion disc structure 

In Table 2 we show several parameters related to the ac- 
cretion structure around the black hole at the end of the 
simulation. The disc masses are computed as before (pa- 
pers II and III), by taking into account the mass which has 
specific angular momentum j > \/&GMtotai /c, so that it will 
remain in orbit around the black hole. This means that only 
a fraction / of the gas mass left in the system at the end of 
the simulation is considered for the accretion tori (usually 
around 70%, the specific value of / for each run is given in 
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Figure 2. continued 



the third column of Table 2). The disc masses are lower at 
lower r (as mentioned above) and at a lower initial mass 
ratio, but there is always at least 0.1 solar masses left in a 
debris torus. Prom the final accretion rate (comparable in 
every case) and the disc mass at t = we infer a rough 
estimate of the lifetime of the disc, which is beteen 30 and 
70 ms. Since it is viscosity that drives the evolution of the 
disc at late times, and hence in this CEise a purely numeri- 
cal effect, these values should be taJcen only as a guideline. 
The true lifetime of the disc depends on the redistribution 
of angular momentum through viscous effects and possibly 
dynamical instabilites. The former would probably make for 



a longer-lived disc, while the latter would tend to act in the 
opposite direction. 

The double ring structure mentioned above gradually 
disappears as the density contrast between the rings drops, 
and the disc becomes more azimuthally symmetric as the 
simulation progresses. By t — tf it is meaningful to take az- 
imuthal averages of quantities such as the density, internal 
energy and specific angular momentum in the disc. These 
are shown in Figure 8 for runs C31 and D31 (the corre- 
sponding plots for the remaining runs arc quite similar). 
The density has a maximum at a characteristic distance ro, 
which is between 50 and 70 km. (Another estimate of the 
lifetime of the disc can be obtained by evaluating Iro/vr-l, 
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Figure 3. Density contours in the orbital plane during the dynamical simulation of the black hole— neutron star binary with initial mass 
ratio q = 0.31 and F = 5/3 (run D31). The orbital rotation is counterclockwise. All contours are logarithmic and equally spaced every 
0.25 dcx. Bold contours arc plotted at logp = —5, —4, —3, —2, —1 (if present) in the units defined in eq. 2. The thick bla<;k arcs bound 
the matter that forms the core (see § 2). The time for each frame is given in the units defined in eq. 1. 



where Vr is the radial (inward) velocity of the fluid at ro. 
The resulting lifetimes Ti„fau are of the same order as those 
shown in Table 2.) The specific internal energies are maxi- 
mum in the inner regions of the discs, and the profiles fiattcn 
out at r « 2ro, at about u/1000 = 2 x 10^'', equivalent to 
3 X 10^® erg g~^, or 3 MeV nucleon"^. The rotation curves 
axe subkeplerian, indicating that pressure support is impor- 
tant. The slight increase in the curves of specific angular 



momentum seen in Figure 8b at r ~ 10 (particularly for 
r = 2) are due to the persistent outer ring in the accre- 
tion disc (see Figure 7). The slow drop in specific angular 
momentum for r > 10 marks the edge of the accretion disc 
eit t = tf. We note that in the inner regions, where the or- 
bits have been circularized through the dynamical process 
mentioned above, the rotation curves axe always close to, 
and below the Keplerian value. The torus does not have a 



© 2001 RAS, MNRAS 000, 1-21 



10 W. H. Lee 






10 15 

Figure 3. continued 




constant distribution of specific angular momentum j, even 
immediately after being formed. 

Since one of our main motivations for this lino of work 
has been the production of cosmological gamma ray bursts 
from theses systems, we show as usual the distribution of 
matter along the rotation axis. This is the region where the 
densities are lowest, and from which a relativistic fireball 
could possibly emerge from the system and produce a GRB. 
This could possibly be powered by neutrino emission from 
the disc and subsequent pair production and expansion, or 
through the Blandford & Znajek (1977) mechanism, by tap- 



ping the energy stored in the spin of the black hole (see 
Table 2). In order for this to occur, the baryon loading must 
be small (on the order of W~^Mq), so that the expansion 
can occur at the required Lorentz factors F > 10^ (Meszaros 
& Rees 1992; Meszaros & Rccs 1993). In Figure 9 we show 
the baryon contamination along the rotation axis by plotting 
the amount of mass enclosed in a cone with opening angle 
above and below the black hole and along the rotation 
axis (see also the last three columns in Table 2). Clearly, 
only modest coUimation is required (of about 10°) to stay 
below the 10~^Mq threshold mentioned above. This is a 
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Figure 4. Separation between the black hole and the centre of mass of the core (see § 2) for (a) F = 2 and (b) F = 5/3 (runs C50 and 
D50 — solid lines; runs C31 and D31 — long— dashed lines; runs C31S and D31S — short- dashed lines; runs C20 and D20-dotted lines). For 
q = 0.31 there are two curves in each frame, corresponding to runs initiated with a spherical polytrope and an irrotational Roche-Riemann 
ellipsoid. In both cases, the one that decays faster corresponds to the former condition. The monotonically decaying curves correspond 
to point-mass binaries with the same initial separation and mass ratio, evolving through gravitational wave emission, computed in the 
quadrupole approximation. 




Figure 5. Mass accretion rate onto the black hole for (a) F = 2 and (b) F = 5/3. The curves are plotted only for t < 100, at later times 
there is little further evolution as M decreases monotonically. 

fact that has become clearer in our simulations as we have rmmber of particles by an order of magnitude, the spatial 
increEised the resolution, from 8,000 SPH particles for the resolution incresises by a factor of ~ 2.15). 
majority of the runs shown in paper I to the 80,000 parti- 
cles used for the simulations shown here (in increasing the 
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Figure 7. Density contour plots at t = tf for runs C31 {a,b) and D31 {c,d) in: (a,c) the orbital plane; (b,d) the meridional plane 
shown by the black line in panels (a,c). All contours are logarithmic and equally spaced every 0.25 dex. Bold contours are plotted at 
logp = —6, —5, —4 (if present) in the units defined in eq. 2. 



Table 2. Accretion disc structure. In the last three columns, 0-n is the half-angle of a cone above the black hole and along the rotation 
axis of the binary that contains a mass M = 10~". The mass is given in units of 1.4 Mq, and time is measured in the units defined in 
equation 1. 



Run 


Mdisc 


/ 


^max 




^ejected 




Jbhc/GM|jj 


e-3 


e-4 


e-5 


C50 


0.181 


0.71 


0.054 


3-10-* 


0.48-10-3 


603 


0.334 


41 


25 


15 


C31 


0.172 


0.71 


0.057 


3-10-* 


10.20-10-3 


573 


0.234 


46 


30 


18 


C31S 


0.179 


0.73 


0.053 


3-10-4 


11.51-10-3 


596 


0.232 


48 


30 


18 


C20 


0.138 


0.59 


0.065 


3-10-4 


6.97-10-3 


460 


0.162 


52 


38 


25 


D50 


0.159 


0.19 


0.067 


310-4 


0.21-10-4 


530 


0.343 


35 


19 


12 


D31 


0.141 


0.63 


0.074 


310-4 


0.80-10-4 


470 


0.242 


43 


26 


12 


D31S 


0.144 


0.64 


0.072 


3-10-4 


2.14-10-4 


480 


0.241 


46 


26 


13 


D20 


0.086 


0.50 


0.094 


3-10-4 


0.05-10-4 


286 


0.172 


51 


35 


22 



© 2001 RAS, MNRAS 000, 1-21 



Black hole-neutron star coalescence 13 



















//w 






— r=2 - 




— r=5/3 







20 40 60 80 100 

t 



Figure 6. Energies in the system as a funcion of time for q = 
0.31 (runs C31 and D31). The kinetic (T), total internal (U), 
gravitational potential (W) and total (E) energies are given in 
units of 3.8 x 10^^ 

4.3 Ejected mass 

During the initial encounter, a tidal tail of material stripped 
from the neutron star is formed in every dynamical run. 
This tail has been observed before for a stiff equation of 
state (paper III), where it persists as a well defined struc- 
ture throughout the simulations. For the eases studied here, 
it is present for F = 2 (see Figure 10), but essentially dis- 
appears as a coherent structure at late times for F = 5/3, 
as the density drops and the tail expands. We have cal- 
culated the amount of dynamically ejected mass for every 
run as before, computing the total mechanical energy (ki- 
netic -l-gravitational potential) of the fluid, and counting as 
ejected that mass for which it has a positive value aX t = tf. 
There are two distinct categories of ejected mass during the 
simulation. The first (type I) corresponds to matter dynam- 
ically ejected from the system, and can be found in the or- 
bital plane, at the tips of the tidal tails formed during the 
disruption of the star at early times {t < 30 — 40, see pan- 
els (b)-(d) in Figures 2 and 3, and Figure 10). The second 
(type II) comes from the surface of the accretion disc, and 
is ejected from the system at later times {t > 70). Ejected 
matter of type II only appears in a significant amount for 
the runs with F — 5/3, and is mainly of numerical origin, 
as it is the equation of state that models more compressible 
gas. This means that the fluid expands to occupy a larger 
volume than for less compressible equations of state. Thus 
for a given number of SPH particles, the spatial resolution is 
lower (i.e. the smoothing length h is larger), particularly at 
the edge of the matter distribution, and the effects of heat- 
ing due to the artificial viscosity can be more pronounced. 
It was not mentioned in paper III simply because no resolv- 
able amount of mass weis ejected in this fashion. For F = 2 
it amounts to only a tiny fraction of the total ejected mass 
(and a handful of particles). For F = 5/3 however, this is 
no longer true. In fact, most of the ejected matter is type II 
in this case. We have not counted it in the values tabulated 



in the sixth column of Table 2, keeping only type I ejected 
matter. Including the internal energy u of the fluid does not 
alter the results given, since the gas coming from the tips 
of the tidal tails has not been subjected to strong compres- 
sion and heating, as it was never part of the accretion torus 
around the black hole. 

For the softer equation of state, mass ejection is strongly 
suppressed, by approximately two orders of magnitude. This 
effect is so strong that for run D20 only about 20 SPH par- 
ticles leave the system. This is similar to what was observed 
by Rosswog et al. (1999) in the case of double neutron star 
mergers. The underlying reason is that as the adiabatic in- 
dex is lowered, the star becomes more centrally condensed, 
and thus the gravitational potential well becomes progres- 
sively deeper. For polytropes, the gravitational potential and 
the density are related by $ = — i<rrF/(F — 1)^"^"^ where 
Kr depends on the value of F. This gives $c(F = 2) = 
0.74<I'c(r — 5/3), at a constant stellar mass and radius. This 
alone makes it more difficult to extract gas from the stellar 
potential well, through the gravitational interaction with the 
black hole during coalescence and eject it from the system. 
In all simulations, we see that the gas that is dynamically 
ejected (type I) comes from the surface layers of the star. 
So if the potential well is deeper, less matter is available for 
this sort of ejection, all other things being equal. There are 
at least three more effects that enhance this result and tend 
to decrease the amount of ejected mass at higher compress- 
ibilities. The first is that more violent events (as measured 
for example by the departure from point-mass behavior at 
small separations, see Figure 4) tend to eject more matter. 
Since these deviations are driven precisely by hydrodynam- 
ical effects, their influence is reduced at low F. Second, as 
pointed out above, the ejected matter comes from the sur- 
face layers of the star, and thus from regions that are at 
lower density at low values of F, making for less total mat- 
ter available for ejection. Third, as can be seen in Figure 1, 
there is less total angular momentum in the system as F is 
decreased (also due to a greater degree of central condensa- 
tion in the star), and so it will be more difficult for matter 
to escape the system in that case. 

The combination of the effects mentioned above makes 
for a dramatic drop in the value of M^jected given in Table 2 
as a function of F. The transition is sharp, due to the gradual 
increase in central condensation of the star, and in particular 
to the qualitative change in the mass-radius relationship 
that occurs at F = 2. 

4.4 Emission of gravitational waves 

The waveforms and luminosities arc calculated in the 
quadrupole approximation from the values of the reduced 
moment of inertia tensor, and its time derivatives, see e.g. 
Finn (1989) and RS92. One polarization of the waveforms 
is shown in Figure 11 for runs C31 and D31, compared with 
the result for point masses decaying in the same approxi- 
mation, and the corresponding luminosities are plotted in 
Figure 12. The results are very similar for the dynamical 
runs with a different initial mass ratio (runs C50, D50, C20 
and D20). 

The tidal disruption of the neutron star in every run 
(irrespective of the mass ratio) after the first episode of 
mass transfer following periastron passage, makes the am- 
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Figure 8. (a) Azimuthally averaged profiles for the density p and internal energy u (w/1000 is plotted) for runs C31 (F = 2) and D31 
(r = 5/3) in the equatorial plane at t = tf. (b) Specific angular momentum j in the equatorial plane for the same runs as in (a). The 
monotonically increasing curve corresponds to that of a Keplerian accretion disc around a black hole of the same mass (the mass of the 
black hole at t = tf for runs C31 and D31 differs by less than 1 per cent). 



A50 




plitudes of the waveforms drop abruptly, and practically to 
zero, as the accretion torus is formed and becomes ever more 
azimuthally symmetric. Upper bounds for the final ampli- 
tude (at t — tf) arc shown in Table 3, where wo show 
the maximum and final amplitudes for the waveforms, the 



peak luminosity and the total energy radiated away by the 
system, and the efficiency of gravitational wave emission 
e = AE / Mtotnic^ ■ For reference, Lmax = 1 (in the units 
given in the table) corresponds to 3.036 x 10 ' ' erg s^^and 
Ai3 = 10 is equivalent to 3.48 x 10^^ erg. The one-armed 
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Figure 10. Density contours in the orbital plane at t = tf for run 
D31. All contours are logarithmic and equally spaced every 0.25 
dex. Bold contours are plotted at log p = —6, —5, —4 in the units 
defined in eq. 2. The thick black line across the tidal tail divides 
the matter that is bound to the black hole from that which is on 
outbound trajectories. 

spiral arms formed during the coalescence (see Figure 10) 
do not contain enough mass to alter the waveforms signifi- 
cantly (Mtaii ~ 0.05). Since the total mass of the system is 
not the same for each run, but increases as the mass ratio 
is decreased, the peak amplitudes in the waveforms (as well 
as the peak luminositites) axe higher as well for lower q (at 
a fixed value of F). At a fixed value of the initial mass ratio, 
however, one can observe the effect of using a different adia- 
batic index clearly. At higher compressibility (i.e. lower F), 
the maximum amplitudes, peak luminosities, the total en- 
ergy release in gravitational waves and the efficiency of this 
emission arc all higher (sec Table 3). The reason for all these 
trends is the same: the higher the compressibility, the more 
centrally condensed the star is. For F = 5/3, Pc/p ~ 5.99, 
while for F = 2, Pc/p = 3.29 (pc is the central density of the 
star and p is its average density). Thus, it resembles a point 
mass to a greater degree in the case with F = 5/3 than if 
F = 2. It is precisely the hydrodynamical effects associated 
with the star not being a point mass that are driving the 
waveforms and luminosities away from the point-mass re- 
sult and making them decay. One can also see in Figure 11 
that for F = 5/3, the waveform takes longer to begin the 
decay, and stays close to the point-mass result for a longer 
time. 



5 INFLUENCE OF INITIAL CONDITIONS ON 
THE DYNAMICAL EVOLUTION OF THE 
SYSTEM 

As for the results we presented in paper III, there are two 
dynamical runs that have used a spherical star as an initial 
condition, instead of an irrotational Roche-Riemann ellip- 
soid. Both have an initial mass ratio q = 0.31, one for F = 2 
(run C31S) and one for F = 5/3 (run D31S). The initial sep- 



10 




Figure 12. Gravitational radiation luminosity for the same runs 
as shown in Figure 11 (solid lines, run C31; dashed lines, run 
D31). The monotonically increasing curves show the correspond- 
ing result for a point— mass binary with the same initial mass 
ratio and separation, decaying in the quadrupole approximation. 
All quantities are given in geometrized units such that G = c = 1 
(Lo = c^/G = 3.64 X 1089 erg s"!). 

aration n is the same as for runs C31 and D31. The initial 
orbital angular velocity n is that for point-mass binaries, 

given that the axis ratios are a^/ai = 02/01 = 1. The pur- 
pose of these runs is to explore the effect of using initial 
conditions that are far from equilibrium for the calculations 
of dynamical coalescence. Since we have already perfomed 
this type of run for a stiff equation of state in paper III, we 
can also gauge how strong the effects are as a function of the 
compressibility. We remind the reader that, even if an irro- 
tational Roche-Riemann ellipsoid is a better approximation 
to the true configuration of the system before coalescence 
than a spherical star, it is not a self-consistent solution, 
since there are no true equilibrium configurations for such a 
system. This is simply because the emission of gravitational 
radiation is always present, and alters the binary separation 
continuously. A tidal lag angle is always present in the bi- 
nary, because the bulge on the surface of the neutron star 
cannot adjust to the changing gravitational potential instan- 
taneously. This angle remains small at large separations, but 
can become quite large (on the order of 10°, see also Lai, 
Rasio & Shapiro 1994) just prior to coalescence, when the 
emission of gravitational waves makes the potential change 
even faster. This aspect of the coalescence process is greatly 
influenced by the stiffness of the equation of state. 

The strongest effect using a spherical star has on the 
dynamical coalescence is due to the response of the star to 
the instantaneous appearance of the gravitational field of the 
black hole at i = 0. A tidal bulge forms, along the line join- 
ing the two centres of mass. The deformed star has a greater 
total energy Wseif + U (the internal energy decreases, see 
Figure 14, but the star is less bound by gravity, and the net 
effect is to increase the energy), which is taken in part from 
the orbital motion, and so the subsequent decay of the orbit 
is faster than for runs C31 and D31 (see Figure 13). The 
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Figure 11. Gravitational radiation waveforms (one polarization is shown) seen at a distance ro away from the system along the rotation 
axis for runs C31 (a), and D31 (b). The dashed lines show the corresponding curves for a point-mass binary with the same initial mass 
ratio and separation, decaying in the quadrupole approximation. All quantities are given in geometrized units such that G = c = 1. 

Table 3. Gravitational radiation. All quantities are given in geometrized units such that G = c = 1, and Lq = c^/G = 3.64 x 10®^ erg s~^ . 



C50 


3.00 


< 0.01 


0.58 


8.95 


4.15-10-3 


C31 


4.25 


< 0.01 


1.13 


15.38 


5.07-10-3 


C31S 


4.27 


< 0.01 


1.15 


13.79 


4.55-10-3 


C20 


5.80 


< 0.01 


2.10 


24.95 


5.7910-'^ 


D50 


3.19 


< 0.01 


0.87 


14.14 


6.5710-'^ 


D31 


4.55 


< 0.01 


1.71 


22.29 


7.35-10--^ 


D31S 


4.49 


< 0.01 


1.58 


21.86 


7.21-10-3 


D20 


6.32 


< 0.01 


3.43 


38.00 


8.82-10-3 



appearance of the tidal bulge also induces radial oscillations 
in the star, which can be seen in the variations of the total 
internal energy U , plotted in Figure 14 at early times (com- 
pare also with the same curves for a stiff equation of state, in 
Figure 14b of paper III). The oscillations are always present, 
but they are somewhat smaller for the runs initiated with 
ellipsoids. 

We shall focus on the results for F = 2 for the follow- 
ing discussion. There arc slight variations if F = 5/3 that 
we will make clear below. As mentioned above, the separa- 
tion initially decreases faster for run C31S than for run C31. 
However, the minimum separation rmin is slighty greater for 
run C31S (see Figure 13), in contreist to the results shown 
in paper III. This is again because of the response of the 
neutron star to mass loss. By expanding and overflowing its 
Roche lobe further after the initial onset of mass loss, the 
encounter develops faster, and the stellar core is pushed out 
to a larger binary separation before approaching the black 
hole any further. This makes the peak accretion rate lower, 
the final disc mass higher, and the Kerr parameter of the 
black hole at t = tf marginally lower (see Table 2). The 
gravitational radiation signal is also affected by the initial 



condition, as can be seen in Table 3. The faster orbital decay 
gives a higher peak amplitude and luminosity (these quan- 
tities depend on the second and third time derivatives of 
the moment of inertia respectively), but a less energetic and 
efficient event, because it is more brief. 

For F = 5/3, the effect on the disc parameters given in 
Table 2 is the same as for F = 2. However, there are qualita- 
tive differences in the way the gravitational radiation signal 
is affected. The energetics and efficiency of the events vary in 
the same way for runs D31 and D31S than for runs C31 and 
C31S, but the trends axe reversed as far as the peak ampli- 
tudes and luminosities are concerned. The reason for this is 
that there are two important factors determining the ampli- 
tude (and hence luminosities) of the gravitational radiation 
waveform: the time derivatives of the moment of inertia, and 
the mass ratio and separation. Inspection of Figure 13b re- 
veals that the differences between runs D31 and D31S are 
small indeed. Essentially, the decay is accelerated by using a 
spherical star, but not nearly fast enough to compensate for 
the fact that the system attains a greater minimum separa- 
tion. Thus the peak amplitudes and luminosities are lower 
for run D31S than for run D31. 
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Figure 13. Separation between the black hole and the centre of mass of the neutron star core as a function of time for (a) runs C31 
and C31S and (b) runs D31 and D31S. The monotonically decaying lines in each frame are the result for a point-mass binary decaying 
through gravitational wave emission, in the quadrupole approximation. 
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Figure 14. Total internal energy U of the neutron star as a 
function of time for runs C31, C31S, D31 and D31S. 



6 SUMMARY, DISCUSSION AND 
CONCLUSIONS 

We have presented the results of three-dimensional dynam- 
ical simulations of the coalescence of a black hole with a 
neutron star, using Smooth Particle Hydrodynamics. The 
black hole is modeled as a Newtonian point meiss with an 
absorbing boundary at the Schwarzschild radius rsch ~ 
2GMbh /(? , and the neutron star is taken as a cold poly- 
trope with adiabatic index F = 2 and T = 5/3. The spatial 
resolution of the results presented here is the highest we 



have used to date, with N ~ 80, 000 SPH particles mod- 
eling the initial neutron star. Dynamical runs with initial 
meiss ratios ranging from q — Mns/Mbh = 0.5 to q = 0.2 
were performed. Given that tidal locking is not expected in 
these systems (Bildsten & Cutler 1992; Kochanek 1992), we 
have used initial conditions that correspond to irrotational 
binaries in equilibrium, approximating the neutron star as 
a compressible tri-axial ellipsoid, following the method of 
LRSb. The dynamical simulations arc begun when the sys- 
tem is on the verge of initiating mass transfer, and followed 
for approximately 23 ms. 

The binary separation decreases as a result of angular 
momentum losses to gravitational radiation, and the neu- 
tron star overflows its Roche lobe within one orbital period 
after the dynamical simulation is started. Irrespective of the 
initial mass ratio and of the value of the adiabatic index, 
this mass transfer episode leads to complete tidal disruption 
of the star on an orbital time-scale. A massive accretion 
disc forms around the black hole, containing a few tenths 
of a solar mass (see Table 2). A single spiral arm appears, 
from material moving through the outer Lagrange point, 
farthest from the black hole. Initially, the accretion torus 
has a complicated structure, with a double ring present (see 
Figures 2 and 3), as the accretion stream collides with it- 
self and circularizes the orbits of the fluid in the disc. As 
the simulation progresses, the disc becomes more and more 
azimuthally symmetric. The peak densities and specific in- 
ternal energies in the discs at the end of the simulations are 
on the order of lO^^g cm~"^ and 10^^ erg g~^ respectively 
(or about 10 MeV nucleon"^). All discs have a low degree of 
baryon contamination along the rotation axis, directly above 
and below the black hole (less than 1O~^M0 are contained 
within approximatlely 10° of the rotation axis). The gravitar 
tional radiation signal reflects the nature of the encounter. 
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with the amplitude of the waveforms dropping practically 
to zero soon after the star is tidally disrupted. Some mass 
{Mejected — 10~^ Mq at most, 866 Table 2), found in the 
outer parts of the tidal tail formed during the initial episode 
of mass transfer, has enough mechanical energy to bo dy- 
namically ejected from the system during coalescence. We 
find that the amount of ejected mass is sensitive to the value 
of the adiabatic index, with a sharp drop (by more than two 
orders of magnitude) ocurring as it decreases below F = 2. 

In paper II we showed the results of dynamical calcu- 
lations of coalescence that used tidally locked binaries with 
an adiabatic index T = 5/3. Thus the effect of using an 
irrotational initial condition can be gauged by directly com- 
paring those results with the present ones. The runs shown 
in paper II also included the effects of gravitational radiation 
reaction in the quadrupole approximation for point masses, 
applying it to the whole star, whereas we have now done it by 
identifying the self-bound core of the neutron star. Qualita- 
tively the outcome of the coalescence is the same for irrota- 
tional and tidally locked systems, but there are quantitative 
differences. These arise because the encounter in the case of 
a tidally locked binary is more gentle, with the separation 
decreasing at a slower rate once hydrodynamical effects be- 
come important. The details can be seen by comparing the 
results given for run D in paper II (tidally locked, T = 5/3, 
initial separation = 3.60 and initial mass ratio q = 0.31) 
with those for run D31 shown here. In the former run, the 
initial peak accretion rate is lower {Mmax = 0.045), the final 
disc mass is higher {Mdisc = 0.226), and the Kerr parameter 
of the black hole is slightly lower (a = 0.222) than for run 
D31 (where Mmax = 0.074, M^sc = 0.141 and a = 0.242), 
all consistent with a less violent encounter after Roche lobe 
overflow. The accretion disc itself is not only more massive, 
but is located at a larger radius, due to the higher value of 
total angular momentum available in synchronized systems. 
This can be soon by locating the maximum in the density 
(which is at r ~ 8 for the tidally locked case, see Figure 10b 
in paper II, and at r ~ 4 for run D31), and the point at 
which the distribution of specific angular momentum j flat- 
tens, marking the outer edge of the disc (at r ~ 13 for run D 
in paper II, and at r ~ 10 for run D31). This also makes the 
maximum densities in the disc greater by at least a factor of 
five in the irrotational case. The tidal tail is almost nonex- 
istent as a large-scale coherent structure for run D31, but 
can be seen clearly in the synchronized case (see Figure 11a 
in paper II). This makes for a lower amount of dynamically 
ejected mass in the irrotational case (by a factor of 200), 
and is due to the lower total angular momentum contained 
in the system, as mentioned above. The exact factor in this 
case remains somewhat uncertain, since it is sensitive to the 
implementation of gravitational radiation reaction, which is 
slightly different in paper II and this work, as mentioned 
above. In paper II we quoted the dynamically ejected mass 
as that which had a positive total energy, including the in- 
ternal energy u. Analysis of those simulations reveals that 
they did contain both type I (cold, dynamically ejected gas) 
and type II ejected matter (see section 4.3). However, as 
mentioned above, there was much more dynamically ejected 
matter than in run D31, and type II matter amounted only 
to « 5 per cent of the total. Finally, the gravitational radi- 
ation signal is affected for the same reasons, with the irro- 
tational case producing a higher maximum amplitude, peak 



luminosity, and total radiated energy (the differences are of 
2, 14 and 7 per cent respectively). 

The present results and those given in paper III allow 
us to observe general trends for all monitored quantities in 
irrotational systems, as the adiabatic index is lowered from 
r = 3tor = 5/3 (see Tables 2 and 3 in paper III and in this 
work). They can be summarized as follows. As the compress- 
ibility increases, the peak accretion rate increases, the disc 
mass drops (here we exclude the results for F = 3 since that 
case did not always imply the complete disruption of the 
neutron star), the black hole has greater spin, the peak am- 
plitude, luminosity and efficiency of gravitational wave emis- 
sion increase, the disc lifetime decreases, and the minimum 
separation attained by the binary before tidal disruption is 
smaller. This last fact implies that the maximum frequency 
emitted by the system in gravitational waves is higher at 
lower F, and can be seen in Figure 15, where we show the 
energy spectrum of the gravitational wave signal for differ- 
ent values of F at a fixed mass ratio q = 0.31. When the 
binary separation is large compared with the stellar radius, 
the spectrum is close to that for a point-mass binary, with 
dE/df oc f~^^^. When the system becomes dynamically un- 
stable, either through tidal effects (for low compressibility) 
or because of runaway mass transfer (for high compressibil- 
ity), the power drops abruptly. This occurs at a characteris- 
tic frequency fdyn, which increases from ~ 700 Hz to 1 kHz 
as the adiabatic index decreases from F = 3 to F = 5/3. 

Each of these consequences can be traced to the degree 
of central condensation of the neutron star, to its mass- 
radius relationship and hence to the way it responds to mass 
loss upon overflowing its Roche lobe. The magnitude of the 
changes in the variables mentioned above is not more than 
a factor of two. The one variable that is greatly affected, 
especially at low values of F, is the total amount of ejected 
mass (see section 4.3). As mentioned in the introduction, the 
mass ejected from this type of system might be a source of 
heavy elements, if the r-process occurs, and could contribute 
significantly to the observed galactic abundances. Our nu- 
merical treatment of the coalescence does not allow us to 
explore nuclear reactions, but merely estimate how much 
matter might leave the system. We refer the reader to the 
work of Rosswog et al. (1999; 2000) and Freiburghaus et 
al. (1999b) for a detailed thermodynamical and nuclear net- 
work calculation, in the case of double neutron stax binares. 
The main point in this respect in our calculations is that (i) 
ejection is greatly suppressed, and practically eliminated, if 
the equation of state is very soft and (ii) irrotational systems 
eject less mass than tidally locked ones, by about one order 
of magnitude. 

The use of accurate equilibrium initial conditions is im- 
portant in dynamical simulations, since an initial perturba- 
tion at the start of the calculation can propagate and affect 
the evolution of the system. Using spherical neutron stars 
for one of our chosen mass ratios, q = 0.31, we have explored 
this effect for irrotational binaries, for values of F ranging 
from 3 to 5/3. We find that the qualitative aspect of the co- 
alescence is unaffected, but that quantitative changes occur, 
all due to the instantaneous appearance of a tidal bulge on 
the surface of the neutron star as the sirrmlation begins. The 
effect of this perturbation is largest at low compressibility, 
since that is when a larger portion of the stellar mass is close 
to the surface, and tidal effects are more pronounced. As one 
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Figure 15. Gravitational waves energy spectrum dE/df for dy- 
namical runs using irrotational binaries with initial mass ratio 
q = 0.31 and F = 3 (solid line, run A31 from paper III), F = 2 
(dashed line, run C31) and F = 5/3 (dotted line, run D31). The 

downward sloping line is the result for a point-mass binary with 
the same mass ratio, with dE/df oc /~^/'^. The increased power 
at / f» 300 Hz for F = 3 corresponds to the return to low fre- 
quencies after the initial mass transfer episode and the survival 
of the binary (see paper III). 



decreases the value of P, the difFcrcnces between runs initi- 
ated with spheres (runs A31S and B31S in paper III and runs 
C31S and D31S in this work) and those that used tri-axial 
ellipsoids (A31, B31, C31 and D31) become less important. 

We have used a polytropic equation of state for our 
study in order to use the compressibility as a free parame- 
ter. Clearly it is an oversimplification as far as thermody- 
namic details are concerned, but it allows one to explore how 
the system responds globally to this variable. As we have 
seen, the emission of gravitational waves and the amount of 
ejected mass are the two aspects that are most affected by 
varying F. One can make the adiabatic index a function of 
the density, and in this way try to model the neutron star 
in a more realistic manner. This approach has been carried 
out by Rosswog et al. (1999; 2000), although they mainly 
used the equation of state of Lattimer & Swesty (1991) for 
their calculations. We have performed tests using this same 
approximation, and have found that it is the value of F at 
high densities that determines the overall evolution of the 
system (as Rosswog et al. did), thus fixing for example the 
qualitative features of the gravitational wave emission and 
the amount of ejected mass. Realistically, it would appear 
that the equation of state for neutron star matter is such 
that the radius is nearly independent of the mass (Prakash 
& Lattimer 2001), and so adopting a polytropic equation of 
state would require using P = 2. 

In a majority of the dynamical simulations we have per- 
formed, we have found that massive accretion discs form, 
with a few tenths of a solar mass. In all cases when this oc- 
curs, the specific angular momentum can be approximated 
by a power law, with j ocr^ . Regardless of the value of the 
adiabatic index, the initial mass ratio, or the initial distribu- 



tion of angular momentum (tidally locked vs. irrotational), 
we find p ~ 0.4 — 0.45. Thus the discs are sub-Keplerian, 
and are far from having a constant distribution of specific 
angular momentum. This is crucial in the context of gamma- 
ray bursts (sec below), because it has been shown that ac- 
cretion discs around black holes can suffer from a runaway 
instability that destroys them on a dynamical time-scale 
(Abramowicz, Calvani & Nobili 1983). Studies over the past 
two decades have shown that a number of effects can ei- 
ther suppress or enhance it. Among these arc (i) the spin 
of the black hole, (ii) the rotation law in the disc, speci- 
fied as j oc r*", (iii) the effects of general relativity and (iv) 
the self- gravity of the disc. Factors (i-high spin) and (ii- 
increasing the value of p) tend to suppress the instability 
(Wilson 1984; Daigne & Mochkovitch 1997; Abramowicz, 
Karas & Lanza 1998; Masuda, Nishida & Eriguchi 1998; Lu 
et al. 2000), while (iii) and (iv) tend to enhance it (Nishida 
et al. 1996; Nishida & Eriguchi 1996; Masuda, Nishida & 
Eriguchi 1998). We note here that all of these studies as- 
sume a softer equation of state than the ones we have used 
(either using a polytrope with P = 4/3 or a realistic equa- 
tion of state for neutron tori). Since our simulations show 
that the Kerr parameter of the black hole is significant, and 
that the power-law index of the distribution of specific an- 
gular momentum is high, it would appear that these discs 
would not suffer from the aforementioned instability, and 
would thus evolve due to angular momentum transport on 
the much longer viscous (rather than dynamical) time-scale. 
However, our simulations are purely Newtonian, and thus it 
is impossible to include the de-stabilizing effects of general 
relativity. The mass of the discs we find is apparently not 
too high (Masuda & Eriguchi 1997), regarding the criterion 
for self-gravity (the mass ratio qdisc = Mdisc/ Mbh between 
the disc and the black hole at the end of the calculations for 
irrotational binaries ranges between 0.02, for run D20, and 
0.09, for run B50 in paper III). 

The accretion discs always have a baryon-frec re- 
gion along the rotation axis, above and below the black 
hole. This region is clear of matter to a degree (less than 
10~^Mq within approximately 10°) that would not hinder 
the production of a relativistic fireball (Mcszaros & Rees 
1992; Meszaros & Rees 1993), thus powering a cosmolog- 
ical gamma ray burst. The binding energy of the tori is 
~ 10^^ erg (see e.g. Figure 6), and the Kerr parameter of the 
black hole is a ~ 0.3 at the end of the calculations, so the 
energy for the burst could come either from neutrino emis- 
sion from the disc or from the spin of the black hole via the 
Blandford & Znajek (1977) mechanism if the magnetic field 
in the torus is strong enough and threads the black hole. 
The maximum extractable energy in this latter case would 
be ~ EBzlO®^ erg, where esz is the MHD efficiency factor. 
In either case, one would expect the disc to survive for a 
time-scale comparable to the duration of the burst, i.e. on 
the order of seconds. This is why the previously mentioned 
result concerning the power law distribution of angular mo- 
mentum and the accompanying dynamical stability of the 
disc is so important. The short time-scales and rapid vari- 
ability involved in a small (« 100 km) accretion disc around 
a black hole make these systems attractive candidates for 
the central engines of short gamma ray bursts, as we found 
in our preliminary studies (Kluzniak & Lee 1998), and have 
now confirmed in the present series of papers for a wide vari- 



© 2001 RAS, MNRAS 000, 1-21 



20 W. H. Lee 



ety of initial conditions, varying the stiffness of the equation 
of state, the initial mass ratio in the binary and the distri- 
bution of angular momentum in the system. 

We note that the mounting observational evidence in 
favor of massive stars being GRB progenitors (Galama & 
Wijers 2001) does not exclude compact mergers as sources, 
simply because all observed afterglows, from which the infer- 
ences about the environment where the bursts occur come, 
correspond to long bursts. If compact mergers do in fact 
produce GRBs, spectacular confirmation about the nature 
of the source could be obtained through the detection of a 
coincident gravitational wave signal, even if the final coales- 
cence waveform is outside the frequency band of detectors 
such as LIGO. One could observe the final minutes of the 
inspiral phase as the orbital frequency increases, leaves the 
LIGO band, and then search for a coincident GRB. 
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APPENDIX A: COMPUTATION OF THE KERR 
PARAMETER OF THE BLACK HOLE 

The Kerr parameter a of the black hole was incorrectly cal- 
culated for the results presented in papers II and III. In 
this appendix we show explicitly the correct derivation of 
Jgg" , and the corrected values for all the runs performed in 
papers II and III. 

When a gas (SPH) particle crosses the accretion bound- 
ary of the black hole, set at the Schwaxzschild radius rsch = 
2GMbh/c^, we update the mass and velocity of the black 
hole so as to ensure conservation of mass and total linear 
momentum, i.e. 

Mbh = Mbh + mi, (Al) 
and 

Mbh^bh = MbhVbh -I- mm, (A2) 

where primed quantities refer to values after the particle has 
been accreted and removed from the simulation. 

The conservation of total angular momentum reads: 

n X mm + fsH X MbhVbh = rsa x M^rVbh + ^bh". (A3) 

where Jg^" is the spin angular momentum gained by the 
black hole because of the accretion. A fraction of the parti- 
cle's angular momentum contributes to the orbital angular 
momentum of the black hole, and the rest to its spin. In 
practice, we found that the latter term dominates, and that 
our error was due mainly to not taking into account the an- 
gular momentum lost to gravitational waves (which is most 
important in the early stages of the simulation, before the 
neutron star has been disrupted). Table Al shows the cor- 
rect values for the Kerr parameter of the black hole for the 
runs shown in paper II (A through E) and for those pre- 
sented in paper III (A50, A31, A31S, A20, B50, B31S, B31 
and B20). 
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Table Al. Kerr parameter of the black hole at the end of the 
dynamical simulations for the runs presented originally in papers 
II and III 



Run 


r 


q 


./bhc/GA/Ijj 


Reference 


A 


5/3 


1.00 


0.448 


paper II 


B 


5/3 


0.80 


0.409 


paper II 


C 


5/3 


0.31 


0.232 


paper II 


D 


5/3 


0.31 


0.222 


paper II 


E 


5/3 


0.10 


0.097 


paper II 


A50 


3.0 


0.50 


0.339 


paper III 


A31 


3.0 


0.31 


0.226 


paper 111 


A31S 


3.0 


0.31 


0.238 


paper 111 


A20 


3.0 


0.20 


0.156 


paper III 


B50 


2.5 


0.50 


0.339 


paper III 


B31 


2.5 


0.31 


0.244 


paper III 


B31S 


2.5 


0.31 


0.247 


paper III 


B20 


2.5 


0.20 


0.167 


paper III 
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